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The problem of the sphtting of a suspension in bifurcating channels dividing into 
two branches of non equal flow rates is addressed. As observed for long, in particular 
in blood flow studies, the volume fraction of particles generally increases in the high 
flow rate branch and decreases in the other one. In the literature, this phenomenon 
is sometimes interpreted as the result of some attraction of the particles towards this 
high flow rate branch. In this paper, we focus on the existence of such an attraction 
through microfluidic experiments and two-dimensional simulations and show clearly that 
such an attraction does not occur but is, on the contrary, directed towards the low flow 
rate branch. Arguments for this attraction are given and a discussion on the sometimes 
misleading arguments found in the literature is proposed. Finally, the enrichment in 
particles in the high flow rate branch is shown to be mainly a consequence of the initial 
distribution in the inlet branch, which shows necessarily some depletion near the walls. 

Key Words: Particle/fluid flow; Microfluidics; Blood flow 



1. Introduction 

When a suspension of particles reaches an asymmetric bifurcation, it is well-known 
that the particle volume fractions in the two daughter branches are not equal; basically, 
for branches of comparable geometrical characteristics, but receiving different flow rates, 
the volume fraction of particles increases in the high flow rate branch. This phenom enon, 
sometimes called the Zweifach-Fung effect (see Svanes fc Zweifach 1968 : FunglllQTSI ). has 



been observed for a long time in the blood circulation. Under standard physiological 
circumstances, a branch receiving typically one fourth of the blood inflow will see its 
hematocrit (volume fraction of red blood cells) drop down to zero, which will have ob- 
vious physiological consequences. The expression 'attraction towards the high flow rate 
branch' is sometimes used in the literature as a synonymous for this phenomenon. Indeed, 
the partitioning not only depends on the interactions between the flow and the particles, 
which are quite complex in such a peculiar geometry, but also on the initial distribution 
of particles. 



Apart the huge number of in- vivo studies on blood flow (see Pries et al. i 1996f) for a 



review), many other papers have been devoted to this effect, either to understand it, or to 
use it in order to design sorting or purification devices. In the latter case, one can play at 
will with the different parameters characterizing the bifurcation (widths of the channels, 
relative angles of the branches) , in order to reach a maximum of efficiency. As proposed in 
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Figure 1. (colour online) (a) The two Y-shaped geometries mainly studied in the literature. 
Here Q\ < Q2 and the dashed line stands for the separating streamline between the flows that 
will eventually enter branches 1 and 2 in the absenc e of particles, (b) The T-bifurcation that is 
studied in this paper and also in lChien et al\ (|1985| ) in order to get rid of geometrical effects as 
much as possible. 



trol this phenomenon (see 


Bugharcllo & Hsiao"l964':'ChiGn et a/.' 1985':'Audet & Olbrichtl 


1987 


; Ditchfield & 01brichtil99&: .Roberts & Olbricht.200a., 2006,;. Yang et a/...2006.; .Barber et al. 


2008 


). In- vitro behavior of red blood cells has also attracted some attention (see Dellimore et al. 


1983 


; Fenton et al. 


198,4 ICarr & Wickham 1990l: Yang et al. 200fi: .Taggi et al\ 20071: 


Zheng et al. I2OO8I: 


Fan et al. 


,2008|). The problem of particle flow through an arrav of 



obst acles, which can be somehow considered as similar has also been studied recently 



(see 



Balvin et al\ 



El-Kareh fc SecombI [2OO0I: iDavis et a/.ll2006l : llnglidboogt [Frechette fc Drazeilliooa 
20091) . 

All the latter papers consider the low- Reynolds- number limit, which is the relevant 
limit for applicative purposes and for the biological systems of interest. Therefore, this 
limit is also considered throughout this paper. 



In most studies as well as in in vivo blood flow studies, which are for historical reasons 
the main sources of data, the main output is the particle volume fraction in the two 
daughter branches as a function of the flow rate ratio between them. Such data can be 
well described by empirical laws that still depend on some ad-hoc parameters but allow 



some rough predictions (see iDellimore et allll983l: iFenton et al. 



1985 



, Pries et al. 
2OIOI) . 



Tgii) . 



which have been exhaustively compared recently (see iGuibert et 

On the other hand, measuring macroscopic data such as volume fractions does not allow 
to identify the relevant parameters and effects involved in this asymmetric partitioning 
phenomenon. 

For a given bifurcation geometry and a given flow rate ratio between the two outlet 
branches, the final distribution of the particles can be straightforwardly derived from two 
data: first, their spatial distribution in the inlet; second, their trajectories in the vicinity 
of the bifurcation, starting from all possible initial positions. If the particles follow their 
underlying unperturbed streamlines (as would a sphere do in a Stokes flow in a straight 
channel), their final distribution can be easily computed, although particles near the apex 
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of the bifurcation require some specific treatment, since they cannot approach it as much 
as their underlying streamhne does. 

The relevant physical question in this problem is thus to identify the hydrodynamic 
phenomenon at the bifurcation that would make flowing objects escape from their un- 
derlying streamlines, which would have as a consequence that a large particle would be 
driven towards one branch while a tiny fluid particle located at the same position would 
go to the other branch. 

In order to focus on this phenomenon, we need to identify more precisely the other 
parameters that influence the partitioning, for a given choice of flow rate ratio between 
the two branches. 



(a) The bifurcation oeometru. \Audet k Olbrichd (|l987h and lRoberts fc Olbrichtl (|2003[ l 
made it clear, for instance, that the partitioning in Y-shaped bifurcations depends strongly 
on the angles between the two branches (see figure [T^). For instance, while the velocity is 
mainly longitudinal, the effective available cross section to enter a perpendicular branch 
is smaller than in the symmetric Y-shaped case. Even in the latter case, the position of 
the apex of the bifurcation relatively to the separation line between the fluids going in 
the two branches might play a role, due to the finite size of the flowing objects. 

(&) The radial distribution in the inlet channel. In an extreme case where all the par- 
ticles are centered in the inlet channel and follow the underlying fluid streamline, they 
all enter in the high flow rate branch; more generally the existence of a particle free 
layer near the walls favours the high flow rate branch, since the depletion in particles 
it entails is relatively more important for the low flow rate branch, which receives fluid 
that occupied less place in the inlet branch. The existence of such a particle free layer 
near the wall has been observed for long in blood circulation, under the name of plasma 
skimming. More generally, it can be due to lateral migration towards the centre, which 
can be of inertial origin (high Reynolds number regime) (see Scho n berg fc HinchI 19891 : 
Asmolovlliool lEloot et alhoo4 iKim fc YoolboOSt IYoo fc Kimll2010l ). or viscous one. In 
such a situation of low Reynolds number flow, while a sphere does not migrate transver- 
sally due to symmetry and linearity i n the Stokes equation, deformable obje cts such as 
vesicles ( closed lipid m embranes) (see 'Coup ler et al. 20081 : Kaoui et al. 2009t ) , red blood 
cells (seelSecomb et al. 2007: Bagchi 2007) that exhibit s imilar dynamics as ve sicles (see 



Abkarian et al. 2007 : Vlahovska et al .I l2009|). drops (see 



et all 2007h or elastic capsules (see 



Secomb et al 



Mortazavi fc T rvggvas onl 12000 
2007t lBagch'ill2007c iRisso et al. 



20061 ) ■ might adopt a shape tha t allows lateral m igrat ion. This migration is due to the 
presence of walls Csee l011al[l997l: lAbkarian e t al. 20o 3: ICallens et alll2008l) as well as to 
the non-constant shear rate (see Kaoui et aLll2008tiDanker et al. 2009f ). Even in the case 
where no migration occurs, the initial distribution is still not homogeneous: since the 
barycentre of particles cannot be closer to the wall than their radius, there is always 
some particle free layer near the walls. This sole effect will favour the high flow rate 
branch. 



(c) Interactions between objects. As illustrated in lDitchfield fc Olbrichtl (|1996I ) or lChesnutt fc Marshall 
(2003), interactions between objects tend to smoothen the asymmetry of the distribution, 
in that the second particle of a couple will tend to go in the other branch as the first one. 
A related issue is the study of trains of drops or bubbles at a bifurcation, that completely 
obstruct the channels and whose passage in the bifurcation greatly modifies the pressure 
dist ribution in its vicinity, and thus influences the behaviour of the follow ing element 
fsee lEngl e< aLl 120051: Ijousse et a^.lliool : ISchindler fc Aidarill2008l: ISessoms et al.\\200^ . 



In spite of the huge literature on this subject, but probably because of the applica- 
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tive purpose of most studies, the relative importance of these different parameters are 
seldom quantitatively discussed, although most authors are fully aware of the different 
phenomena at stake. 

As we want to focus in this paper on the question of cross streamline migration in the 
vicinity of the bifurcation, we will consider rigid spheres, for which no transverse migra- 
tion in the upstream channel is expected, that are in the vanishing concentration limit 
and flow through symmetric bifurcations, that is the symmetric Y-shaped and T-shaped 
bifurcations shown on figure [TJ where the two daughter branches have same cross section 
and are equally distributed relatively to the inlet channel. 

Indeed, this rigid spheres case is already quite unclear in the literature. In the following, 
we first make a short review of some previous studies that consider a geometrically 
symmetric situation and thoroughly re-analyze their results in order to detect whether 
the Zweifach-Fung effect they see is due to initial distribution or to some attraction in 
the vicinity of the bifurcation, which was generally not done (section [2]). 

We then present in sections [3] and |4] our two-dimensional simulations and quasi-two- 
dimensional experiments (in a sense that the movement of the three-dimensional objects 
is planar). We mainly focus on the T-shaped bifurcation, in order to avoid as much as 
possible the geometrical constraint due to the presence of an apex. 

Our main result is that there is some attraction towards the low flow rate branch 
(section |4.1[) . This result is then analyzed and explained through basic fluid mechanics 
arguments, which are compared to the ones previously evoked in the literature. 

In a second time, we discuss which consequences this drift has on the flnal distribution 
in the daughter branches. To do so, we focus on what the particles concentrations at 
the outlets would be in the simplest case, that is particles homogeneously distributed in 
the inlet channel, with the sole (and unavoidable) constraint that they cannot approach 
the walls closer than their radius (denominated as depletion effect in the following, see 
figure [TJb)). This is done through simulations, which allow us to easily control the initial 
distribution in particles (section l4.2[) . Consequences for the potential efficiency of sorting 
or purification devices are discussed. We finally come back, in section [531 to some of the 
previous studies found in the literature with which quantitative comparisons can be done 
in order to check the consistency between them and our results. 

Before discussing the results from the literature and presenting our own data, we shall 
introduce useful common notations (see figure [Dd). 

The half-width of the inlet branch is set as the length scale of the problem. The inlet 
channel divides into two branches of width 2a (the case a = 1 will be mainly considered 
here by default, unless otherwise stated), and spheres of radius i? ^ 1 are considered. 
The flow rate at the inlet is noted Qo, and Qi and Q2 are the flow rates at the upper 
and lower outlets (Qo = Qi + Q2)- In the absence of particles, all the fluid particles 
situated initially above the line y — yf will eventually enter branch 1. This line is called 
the (unperturbed) fluid separating streamline, yo is the initial transverse position of the 
considered particle far before it reaches the bifurcation (|?/o| ^ 1 — R). Ni and N2 are 
the numbers of particles entering branches 1 and 2 by unit time, while Nq = Ni + N2 
have entered the inlet channel. The volume fractions in the branches are = VNi/Qi, 
where V is the volume of a particle. 

With these notations, we can reformulate our question: if j/o = Uf, does the particle 
experience a net force in the y direction (e. g. a pressure difference) that would push 
it towards one of the branches, while a fluid particle would remain on the separating 
streamline (by definition of y/)? If so, for which position ?/q does this force vanish, so 
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that the particle follows the streamlines and eventually hits the opposite wall and reaches 
an (unstable) equilibrium position? If Qi ^ Q2 and j/g < j//, then one will talk about 
attraction towards the low flow rate branch. 
Following these notations, we have: 

TVi = / n{y)ul{y)dy, (1.1) 

•'Vo 

Qi^ ( u^{y)dy, (1.2) 

where n(jj) is the mean density in particles at height y in inlet branch, u* and Ux are 
respectively the particles and flow longitudinal upstream velocities. iVo and Qq are given 
by the same formula with j/g = y/ = —1. 

The Zweifach-Fung effect can then be written as follows: if Qi/Qo < 1/2 (branch 1 
receives less flow than branch 2) then Ni/Nq < Qi/Qo (branch 1 receives even less par- 
ticles than fluid) or equivalently $1 < (f>g (the particle concentration is decreased in the 
low flow rate branch). 



2. Previous results in the literature 

In the literature, the most common symmetric case that is considered is the Y-shaped 
bifurcation with daughter branches leaving the bifurcation with a 45° angle relatively to 
the inlet channel, and cross sections identical as the one of the inlet channel (figure H^) 
(see 'Audet & 01bricht"l987'; 'Pitchf ield fc 01bricht|[l99l [Roberts fc Olbrichtll200d l200a 
[Vang et al. 2006; Barber et al. 20 j) . The T-shaped bifurc ation (figurelDa) has attracted 
little attention fsee lYen fc Fungi Il978t IChien et aLlll985h . AU studies but lYen fc Fungi 
( 19781 ) show res ults for rigi d spherical parti cles, while some results for deformable parti- 
cles are given in Yen fc Fun g (1978) and B arber et al. ( 20081 ). Explicit data on a possible 



attraction towards one branch are scarce as fhey can only be found in a recent two- 
dimensional simulations paper (see iBarber et al" 2008*) . In three other papers, dealing 
with two-dimensional simu lations (see I Audet & : Olbricht 1987) or expe riments in square 



cross section channels (see Roberts fc Olbrichtlbood [Yang a/.ll2006l ). the output data 
are the concentrations at the outlets. In this section, we re-analyze their data in order 
to discuss the possibility of an attraction toward s one branch. Ex periments in circu- 
lar cross section charinels were also developed (see Yen fc Fund Il97i 



Chien et al. 1985 



Pitchfield fc 01brichtl[l995 iRoberts fc Olbrichtll2003n , on which we comment in a second 



time. 

In the two-dimensional simulations presented in lAudet fc Olbrichti (|l987n . some trajec- 
tories around the bifurcation are shown, however the authors focused on an asymmetric 
Y-shaped bifurcation. In addition, some data for Ni/Nq in a symmetric Y-shaped bifur- 
cation and R = 0.5 are presented. Yang et al. ran experiments with balls of similar size 
(R = 0.46) in a symmetric Y-shaped bifurcati on with s quare cross section and also showed 
data for Ni/Nq as a function of Qi/Qo (see lYang et al. 2006). Exp erime nts with l arger 
balls ( R = 0.8) in square cross section channels were carried out in iRoberts fc Olbricht 
( 20061 ). Once again, the output data are the ratios Ni/Nq. In both experiments, the 
authors made the assumption that the initial ball distribution is homogeneous, as con- 
sidered also in the simulation paper by Audet and Olbricht. In all the latter papers, 
although the authors are sometimes conscious that the depletion and attraction effects 
might screen each other, the relative weight of each phenomenon is not really discussed. 
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Figure 2. Comparison between data from literature and theoretical distribution under the 
assumption of no attraction, which would indicate some previously unseen attraction towards 
the low flow rate branch. Rigid spheres distribution Ni/Nq is shown as a function of flow 
distribution Qi/Qo in a symmetric Y-shaped bifurcation and an homogeneous distribution at 
the inlet (but the una voidable depletion effect). Symbols: d ata extracted from previous papers: 
(■ lYang et al\ ||2006D . figure 3, experimen ts, R = 0.46. (AlAudet fc Olbrichtl ifToST), figure 8, 
two-dimensional simulations, R — 0.5. (□) iRoberts fc Olbrichtl (j2006l ). figure 5A. experiments, 
R = 0.8. Dotted and full lines: theoretical distribution for R = 0.48 and R — 0.8 in case the 
particles follow their underlying streamline (j/o = yf- no-attraction assumption) and u* given 
by our simulations. Dashed line: fluid distribution {Ni/No = Qi/Qo)- 



However, Yang et al. consider explicitly that there must be some attraction towards the 
high flow rate branch an d give some qualitative argum e nts for it. T his opinion, initially 
introduced by Fu ng (see Fung 1973t Yen fc Fung 1978 : Fuii3 1993), is widely spread in 



the literature (see IeI Kareh fc Secombll2000Hjaggi et al\\200l[ iKersaudv-Kerhoas et al 



20ld) . We shall come back to the underlying arguments in the following. 



In f igure m we present the data of Ni/Nq as a fun ction of Qi / Qo taken from Audet fc Olbricht 
( 1987t ) for R = 0.5 (two-dimensional sim ulations), lYang et all (^2006) for R = 0.46 (exper- 
iments) and lRoberts fc Olbrichtl ( 20061 ) for R — 0.8 (experiments). It is very instructive 
to compare these data with the corresponding values calculated with a very simple model 
based on the assumption that no particular effect occurs at the bifurcation, that is, the 
particles follow their underlying streamline (no- attraction assumption). To do so, we 
consider the two-dimensional case of flowing spheres and calculate the corresponding Ni 
according to equation (jl.ip . The no- attraction assumption implies that Vq — Vf and, 
as in the considered papers, the density n{y) is considered constant for |y| ^ 1 — R. 
The particles velocity u* is given by our simulations presented in section 14.21 Since we 
consider only flow ratios, this two-dimensional approach is a good enough approximation 
to discuss the results of the three-dimension experiments, as the fluid separating plane 
is orthogonal to the plane where the channels lie; moreover, the position of this plane 
differs only by a few percent from the position of the separating line in two dimensions. 
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In all curves, it is seen that, if Qi/Qo < 1/2, then Ni/Nq < Qi/Qq, which is precisely 
the Zweifach-Fung effect. Note that this effect is present even under the no-attraction 
assumption: as already discussed, the sole depletion effect is sufficient to favour the high 
flow rate branch. 

Let us first consider spheres of medium size {R ~ 0.48: Audet and Olbricht / Yang et 
ai). If we compare the data from the literature with the theoretical curve found under 
the no-attraction assumption, we see that the enrichment in particles in the high fiow 
rate branch is less pronounced in the simulations by Audet and Olbricht and of the 
same order in the experiments by Yang et al. Therefore, we can assume that in the 
two-dimensional simulations by Audet and Olbricht, there is an attraction towards the 
low flow rate branch, which lowers the enrichment of the high flow rate branch. The case 
of the experiments is less clear: it seems that no peculiar effect takes place. 

The R = 0.8 case is even more striking: under the no-attraction assumption, we can 
see that for Qi/Qo < 0.35, A^i = because yf > 1 — R and no sphere can enter the 
low flow rate branch. In the meantime, a non negligible amount of particles are found 
to enter branch 1 for Qi/Qo < 0.35 by Roberts and Olbricht in their experiments (see 
figure [2]). It is clear from this that there must be some attraction towards the low flow 
rate branch. 



For channels with circular cross sections, the data found in the lite rature do no t all te ll 
the same story, although spheres of similar sizes are considered. In lChien et ali (|l985r ). 
R = 0.79 spheres are considered in a T-shaped bifurcation. The Y-shaped bifurcation 
was considered twice by the same research grou p, with very similar sphere s: R = 0.8 (see 
Ditchfield fc 01brichtlll996l ) and R = 0.77 (see lRoberts fc Olbricht! l2003[ ). In a circular 



cross section channel, the plane orthogonal to the plane where the channels lie, parallel 
to the streamlines in the inlet channel and located at distance 0.78 from the inlet channel 
wall corresponds to the flow separating plane for Qi/Qo — 0. 32. At low c o ncent rations, 
very few spheres are observed in branc h 1 for Qi/Qq < 0.32 in lChien et al. i 19851 ) (figure 
3D) and Ditchfield fc Olbrichtj (Il996h (figure 3), in agreement with a no-attraction as- 
sumption. In Chien et'oll (1985), the authors also show their data can be well described 
by the theoretical curve calculated by assuming the particles follow their underlying 
streamlines. In marked contrast with these results, a considerable amount of sp heres is 
still observed in branch 1 in the same situation in [Roberts fc OlbrichtJ (|2003l ^ (fi gure 
4). Similarly, in Ditchfield fc Olbricht ( 19961 ) (figure 4), many particles with R = 0.6 are 
found to enter the low flow rate branch 1 even when Qi/Qo < 0.19, which would indicate 
some attraction towards the low flow rate branch. Thus, in a channel with circular cross 
sectio n, the results are contradictory. In the pioneering work presented in Yen fc Fund 
( 1978h . a T-shaped bifurcation is also considered, with flexible disks mimicking red blood 
cells, but the deformability of these objects and the noise in the data do not allow us to 
make any reasonable discussion. 



More recently, iBarber et al. (|2008h presented simulations of two-dimensional spheres 
with R ^ 0.67 and two-dimensional deformable objects mimicking red blood cells in a 
symmetric Y-shaped bifurcation. The values of as a function of the flow rate ratios and 
the spheres radius are clearly discussed. For spheres, it is shown that j/g < ?// if Qi < Q2, 
that is, there is an attraction towards the low flow rate branch, which increases with R. 
Deformable particles are also considered. However, it is not possible to discuss from their 
data (as, probably, from any other data) whether the cross streamline migration at the 
bifurcation is more important in this case or not: for deformable particles, transverse mi- 
gration towards the centre occurs, due to the presence of walls and of non homogeneous 
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shear rates. This migration will probably screen the attraction effect, at least partly, and 
it seems difficult to quantify the relative contribution of both eff ects. In particular, ?/n 
depen ds on the (arbitrary) initial distance from the bifurcation. In lChesnutt &: Marshall 
( 2009( ). attraction towards the low flow rate branch is also quickly evoked, but considered 
as negligible since the authors mainly focus on large channels and interacting particles. 



Finally, from our new analysis of previous results of the literature (and despite some 
discrepancies) it appears that there should be some attraction towards the low flow rate 
branch, although the final result is an enrichment of the high flow rate branch due to 
the depletion effect in the inlet channel. This effect was seen by Barber et al. in their 
si mulations. On the other han d, if one considers the flow around an obstacle, as simulated 
in El-Kareh fc Secomb ( 2000f) . it seems that spherical particles are attracted towards the 



high flow rate side. 



From this we conclude that the different effects occurring at the bifurcation level are 
neither well identifled nor explained. Moreover, to date, no direct experimental proof of 
any attraction phenomenon exists. In section 14. 1[ we show experimentally that attrac- 
tion towards the low flow rate branch takes place and conflrm this through numerical 
simulations. 

It is then necessary to discuss whether this attraction has important consequences 
on the final distributions i n par ticles in the two daughter channels. This was not done 
explicitly in iBarber et al\ (I2OO8"). It is done in section IT2l where we discuss the relative 
weight of the attraction towards the low flow rate branch and the depletion effect, which 
have opposite consequences, by using our simulations. 



3. Method 

3.1. Experimental setup 

We studied the behaviour of hard balls as a flrst reference system. Since the potential 
migration across streamlines is linked to the way the fluid acts on the particles, we also 
studied spherical fluid vesicles. They are closed lipid membranes enclosing a Newtonian 
fluid. The lipids that we used are in liquid phase at room temperature, so that the 
membrane is a two-dimensional fluid. In particular, it is incompressible (so that spherical 
vesicles will remain spherical even under stress, unlike drops), but it is easily sheared: 
it entails that a torque exerted by the fluid on the surface of the particle can imply a 
different response whether it is a solid ball or a vesicle. Moreover, since vesicle suspensions 
are polydisperse, it is a convenient way to vary the radius R of the studied object. 

The experimental setup is a standard microfluidic chip made of polydimethylsiloxane 
bonded on a glass plate (figure |3]). We wish to observe what happens to an object located 
around position yf that is, in which branch it goes at the bifurcation. In order to de- 
termine the corresponding y^, we need to scan different initial positions around yf. One 
solution would be to let a suspension flow and hope that some of the particles will be 
close enough to the region of interest. In the meantime, as we shall see, the cross stream- 
line effect is weak and requires precise measurement, and noticeable effects appear only 
at high radius i?, typically R > 0.5. With such objects, clogging is unavoidable, which 
would modify the flow rates ratio, and if a very dilute suspension is used, it is likely that 
the region of interest will only partly be scanned. 

Therefore we designed a microfluidic system allowing to use only one particle, that 
would go through the bifurcation with a controlled initial position would be taken 
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Figure 3. (colour online) Scheme of the microfluidic device. The photograph shows the trajec- 
tory of a particle from branch a to branch 1 after having been focused on a given streamline 
thanks to flows from lateral branches b and c. 



back, its position j/o modified, would flow again in the bifurcation, and so on. Moreover, 
we allowed continuous modification of the flow rate ratio between the two daughter 
branches. The core of the chip is the five branch crossroad shown in inset on figure [H 
These five branches have different lengths and are linked to reservoirs placed at different 
heights, in order to induce a flow by hydrostatic pressure gradient. A focusing device 
(branches a,b and c) is placed before the bifurcation of interest (branches 1 and 2), in 
order to control the lateral position of the particle. Particles are initially located in the 
central branch a, where the flow is weak and the incoming particles are pinched between 
the two lateral flows. In order to modify the position i/q of the particle, the relative 
heights of the reservoir linked to the lateral branches are modified. The total flow rate 
and the flow rate ratios between the two daughter branches after the bifurcation are 
controlled by varying the heights of the two outlet reservoirs. Note the flow rates ratio 
also depends on the heights of the reservoirs linked to inlet branches a, b and c. Since the 
two latter must be continuously modified to vary the position yo of the incoming particle 
in order to find ?/q for a given flow rate ratio, it is convenient to place them on a pulley 
so that their mean height is always constant (the resistances of branches b and c being 
equal). If the total flow rate is a relevant parameter (which is not be the case here since 
we consider only Stokes flow of particles that do not deform), one can do the same with 
the two outlet reservoirs. In such a situation, if reservoir of branch a is placed at height 
0, reservoirs of branches b and c at heights ±/io, and reservoirs of branches I and 2 at 
height —H + h and —H — h, the flow rate ratio is governed by setting (h, H) and /iq can 
be modified independently in order to control j/o- Once the particle has gone through 
the bifurcation, height H and the height of reservoir a are modified so that the particle 
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Figure 4. Photographs showing the difTerent positions of a vesicle of radius R — 0.60 starting 
just above and just below its separating line. No clear difference between these two starting 
positions can be seen with eyes, which illustrates the accuracy we get in the measurement of yg . 
Qi/Qo is set to 0.28. 



comes back to branch a, and ho is modified in order to get closer and closer to position 
Vo- Qi/Qo (or, equivalently, y/) is a fmiction of H, h, and the flow resistances of the 
five branches of rectangular cross sections, which are known functions of their lengths, 
widths and thicknesses (seciWhite 1991). The accuracy of the calculation of this function 
was checked by measuring Uq for small particles, that must be equal to y/. 

Note that the length of the channel is much more important than the size of a single 
flowing particle, so that we can neglect the contribution of the latter in the resistance to 
the flow: hence, even though we control the pressures, we can consider that we work at 
fixed flow rates. 

Finally, as it can be seen on figure |4l our device allows us to scan very precisely the 
area of interest around the sought j/q , so that the uncertainty associated to it is very low. 

At the bifurcation level, channels widths are all equal to 57 ± 0.2/im. Their thickness 
is 81 ± 0.3/j,m. We used polystyrene balls of maximum radius 40.5 ± 0.3 /xm in soapy 
water (therefore R ^ 0.71) and fluid vesicles of size R ^ 0.60. Vesicles membrane is a 
dioleoylphosphatidylcholine lipid bilayer enclosing an inner solution of sugar (sucrose or 
gluc ose) in water. Vesicle s are produced following the standard electroformation method 
(see Angelova et aL|[l992l ). Maximum flow velocity at the bifurcation level was around 1 



mm.s ^, so that the Reynolds number Re ~ 10 ^. 

3.2. The numerical model 

In the simulations, we focus on the two-dimensional problem (invariance along the z 
axis). Our problem is a simple fluid/structure interaction one and can be modeled by 
Navier-Stokes equations for the fluid flow and Newton-Euler equations for the sphere. 
These two problems can be coupled in a simple manner: 

• The action of fluid on the sphere is modeled by the hydrodynamic force and torque 
acting on its surface. They are used as the right hand sides of Newton-Euler equations. 

• The action of the sphere on fluid can be modeled by a non-slip boundary conditions 
on the sphere (in the Navier-Stokes equations) . 

However, this explicit coupling can be unstable numerically and its resolution often re- 
quires very small time steps. In addition, as we have chosen to use finite element method 
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FEM (for accuracy reasons) and since the position of the sphere evolves in time we have 
to remesh the computational domain at each time step or in best cases at each few time 
steps. 

For all these reasons we chose another strategy to model our problem. Instead of using 
Newton-Euler equations for modeling the sphere motion and Navier-Stokes equations for 
the fluid flow, we use only the Stokes equations in the whole domain of the bifurcation 
(including the interior of the sphere). The use of Stokes equations is justified by the 
small Reynolds number in our case and the presence of the sphere is rendered by a 
second fluid with a 'huge' viscosity on which we impose a rigid body constraint. This 
type of strategy is widely used in the litera ture with diff erent names e.g. the so called 



FPD (Fluid Particle Dynamics) method (see lTanaka fc Araki 2000 : Pey la .2QQ7) ) but we 



can group them under the generic name of penalty-like methods. The one t hat we use is 
mainly developed by Lefebvre et al. (see Janela et al. 2005^ Lefebvre"2007^ and we can 



find a mathematical analysis of these types of methods in [Maury (2009) . 

In what follows we describe briefly the basic ingredients of the finite element method 
and penalty technique applied to our problem. 

The fluid flow is governed by Stokes equations that can be written as follow: 

- :/Au + Vj3 = in (3.1) 

V-u = 0in%, (3.2) 

u = f on dilf. (3.3) 

Where: 

m u and p are respectively the viscosity, the velocity and the pressure flelds of the 
fluid, 

• rjy is the domain occupied by the fluid. Typically flf — fl \ B if we denote by Q the 
whole bifurcation and by B the rigid particle, 

• 9ri/ is the border of SI/, 

• f is some given function for the boundary conditions. 

It is known that under some reasonable assumptions the problem (I3.1l) - (l3.2p - ()3.3p has a 

unique solution (u,p) £ H^{Qf)'^ x L§(ri/) (see iGirault fc Raviartiil986) . In the sequel 
we will use the following functional spaces: 

L^{n) ^{f:n~^R;[ \f\^ < +c»}, (3.4) 
Llin) = {f e L^{n); [ / = 0}, (3.5) 

H\n)^{f eL^{ny,vf £L\n)}, (3.6) 

Ho{n) = {/ e H^(n):f = on dn}. (3.7) 

As we will use FEM for the numerical resolution of problem p.l|) - p.2p - p.3|) . we need 
to rewrite it in a variational form (an equivalent formulation of the initial problem) . For 
sake of simplicity, we start by writing it in a standard way (fluid without sphere), then 
we modify it using penalty technique to take into account the presence of the particle. In 
what follow we describe briefly these two methods, the standard variational formulation 
for the Stokes problem and the penalty technique. 

3.2.1. Variational formulation 

Let us first recall the deformation tensor r which will be useful in the sequel 

T(u) = i(Vu+(Vu)*). (3.8) 
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Thanks to incompressibility constraint V • u = we have 



Au = 2V-r(u). (3.9) 

Hence, the problem (I3.ip - p.2p - (l3.3l) can be rewritten as follows: find (u,p) S i?^(57/)^ x 
^□(fi/) such that: 

-2i^V-r(u) + Vp = in 17/, (3.10) 

V-u = Oinll/, (3.11) 

u = fonai7/. (3.12) 

By simple calculations (see appendix for details) we show that problem p.lOp - p.lip - 
(|3.12p is equivalent to this one: find (u,p) G H'^{VlfY x -^o(^/) ^^ch that: 

21//" T(u):r(v)-/" pV • v = 0, Vv e i/i(%)2, (3.13) 

I qV •U = 0,WqeLl{nf), (3.14) 

u = fonarj/, (3.15) 

where : denotes the double contraction. 

3.2.2. Penalty method 

We chose to u se the penalty strat e gy in the framew ork of FEM that we will describe 



briefly here (see lJanela et al\ (|2005[ ): iLefebvrd (|2007l ) for more details) 



The first step consists in rewriting the variational formulation (I3.13p - (l3.14p - (l3.15|) by 
replacing the integrals over the real domain occupied by the fluid {ilf — i^\B) by those 
over the whole domain J7 (including the sphere B). Which means that we extend the 
solution {u,p) to the whole domain ft. More precisely, by the penalty method we replace 
the particle by an artificial fiuid with huge viscosity. This is made possible by imposing 
a rigid body motion constraint on the fiuid that replaces the sphere (t(u) = in B). 
Obviously, the divergence free constraint is also insured in B. 

The problem (|3T3l) - (i3T4l) - (j3T4| is then modified as follows: find {u,p) e H^{n)^xLl{n) 
such that: 

2iy f r(u) : t(v) + - / r(u) : r(v) 

- [ • V = o,Vv e i^ol(^7)^ (3.16) 

Jn 

f qV ■U = 0,yqeLl{n), (3.17) 
Jn 

u = f on dil. (3.18) 

Where e 1 is a given penalty parameter. 

Finally, if we denote the time discretization parameter by t„ = nSt, the velocity and 
the pressure at time tn by (u„,p„), the velocity of the sphere at time t„ by V„ and its 
centre position by X„, we can write our algorithm as: 

V„ = ^. ■ ^ / u„ (3.19) 
V olume(B) jg 

X„+i = X„ + 6tYn (3.20) 
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(u„+i,p„+i) solves: 

2!^ f t(u„+i) : t(v) + - / t(u„+i) ; r(v) 
Jn £ JB 

~ [ p„+iV-v = 0,Vvei7i(l])2, (3.21) 
Jn 

/ (7V-u„+i = 0,VgeLg(r!), (3.22) 
Jn 

u„+i = f on dn. (3.23) 
The implementation of algorithm (IXTOt-dX^- d^T^-dHT^-dH:^ is d one by using a 



user- friendly finite element software: Freef em++ (see Hecht fc Pironneau 201Cll ). 

Finally, we consider the bifurcation geometry shown in figure [IJb) and impose no-slip 
boundary conditions on all walls and we prescribe parabolic velocity profiles at the inlets 
and outlets such that, for a given choice of flow rate ratio, Qo = Qi + Q2- For a given 
initial position j/o of the sphere of given radius R at the outlet, the full trajectory is 
calculated until it definitely enters one of the daughter branches. A dichotomy algorithm 
is used to determine the key position j/q. Spheres of radius R up to 0.8 are considered. 

Remark 1. In practice, the penalty technique may deteriorate the preconditionning 
of our underlying linear system. To overcome this problem, one can regularize equation 
I13.22\) by replacing it with this one: 

-eo Pn+iq+ / qV -Un+i ^0,yq e Llifl), (3.24) 
Jn Jn 

where Eq <^ 1 is a given parameter. 



4. Results and discussion 

4.1. The cross streamline migration 

4.1.1. The particle separating streamlines 

In figure [5] we show the position of the particle separating line j/q relatively to the 
position of the fluid separating line y/ when branch 1 receives less fluid than branch 2 
(see flgure[T]D), which is the main result of this paper. For all particles considered, in 
the simulations or in the experiments, we flnd that the particle separating line lies below 
the fluid separating line, the upper branch being the low flow rate branch. These results 
clearly indicate an attraction towards the low flow rate branch: while a fluid element 
located below the fluid separating streamline will enter into the high flow rate branch, a 
solid particle can cross this streamline and enter into the low flow rate branch, providing 
it is not too far initially. It is also clear that the attraction increases with the sphere 
radius R. 

In particular, in the experiments (figure[S^), particles of radius R < 0.3 behave like fluid 
particles. R = 0.52 balls show a slight attraction towards the low flow rate branch, while 
the effect is more marked for big balls of radius R = 0.71. Vesicles show comparable 
trend and it seems from our data that solid particles or vesicles with fluid membrane 
behave similarly in the vicinity of the bifurcation. 

In the simulations (flgure[5lD) we see clearly that for a given R, the discrepancy between 
the fluid and particle behaviour increases when Qi/Qo decreases. On the contrary, in 
the quasi-two-dimensional case of the experiments, the difference between the flow and 
the particle streamlines seems to be rather constant in a wide range of Qi/Qo values. 
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Figure 5. Position of the particle separating line y^. The T-bifurcation with branches of equal 
widths is considered. Branch 1 receives flow from high y values, so j/J < j// for Qi/Qo < 1/2 
indicates attraction towards the low flow rate branch (see also figure [TJj). (a) Data from quasi- 1- 
wo-dimensional experiments and comparison with two-dimensional case for one particle size. 
The two-dimensional and three-dimensional fluid separating lines are shown to illustrate the 
low discrepancy between the two cases, as requested to validate our new analysis of the lit- 
erature in section [2l The horizontal dotted line shows the maximum position yo — 1 — R for 
R — 0.71 spheres. Its intersection with the curve yo{Qi/Qo) yields the critical flow rate ratio 
Qi/Qo below which no particle enters branch 1, the low flow rate branch. This expected crit- 
ical flow rates for the two- and three-dimensional cases are shown by arrows, (b) Data from 
two-dimensional simulations. 
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Finally, for small enough values of Qi/Qo, the attraction effect is more pronounced in 
the two-dimensional case than in the quasi- two-dimensional one, as shown on figure [5][a) 
for R = 0.71. This was to be expected, since this effect has something to do with the 
non zero size of the particle and the real particle to channel size ratio is lower in the 
experiments for a given R, due to the third dimension. In all cases, below a given value 
of Qi/Qo, the critical position j/q would enter the depletion zone i/q > 1 ^ R, so that no 
particle will eventually enter the low flow rate branch. The corresponding critical Qi/Qo 
is much lower in the two-dimensional case than in the experimental quasi-two-dimensional 
situation (see figure [5^). 

4.1.2. Discussion 

T he first argument for some attraction towar ds one branch was initially given by Fung 
(sec Fung 'l973'; 'Yen fc Fund Il978l : iFungl [1993I ) and strengthened by recent simulations 



(see Yang fc Zahn , 20041): a sphere in the middle of the bifurcation is considered (j/o = 0) 
and it is argued that it should go to the high fiow rate branch since the pressure drop 
Pq — P2 is higher than Pq — Pi because Q2 > Qi (see figure [Ijb) for notations). This is 
true (we also found j/q > when Qi < Q2) but this is not the point to be discussed: if one 
wishes to discuss the increase in volume fraction in branch 2, therefore to compare the 
particles and fiuid fluxes N2 and Q2, one needs to focus on particles in the vicinity of the 
fluid separating streamline (to see whether or not they behave like the fiuid) and not in 
the vicinity of the middle of the channel. On the other hand, this incorrectly formulated 
argument by Fung has led to the idea that there must be some attraction towards the 



high fiow rate branch in the vicinity of the fiuid separating streamline (see Yang et 



2006 ) , which appears now in the li terature as a well established fact (see Ijaggi et 



2007t ^Kersaudv-Kerh oas et aZ.|[2010l ). 



In iBarber e t al. (2 008h , Fung's argument is rejected, although it is not explained why. 
Arguments for attraction towards the low fiow rate branch (that is, P2 > Pi on figure 
[ijb)) are given, considering particles in the vicinity of the fiuid separating streamline. The 
authors' main idea is, first, that some pressure difference Pq — Pi builds up on each side 
of the particle because it goes more slowly than the fluid. Then, as the particle intercepts 
a relatively more important area in the low flow rate branch region (yf < y < 1) than 
in the high fiow rate region, they consider that the pressure drop is more important in 
the low flow rate region, so that P2 > Pi- The authors call this effect 'daughter vessel 
obstruction'. 

Indeed, it is not clear in this paper where the particles must be for this argument to 
be valid: at the entrance of the bifurcation, in the middle of it, or close to the opposite 
wall as we could think since their arguments are used to explain what happens in case of 
daughter branches of different widths. Indeed, we shall see that the effects can be quite 
different according to this position and, furthermore, the notion of 'relatively larger part 
intercepted' is not the key phenomenon to understand the final attraction towards the 
low flow rate branch, even though it clearly contributes to it. 

To understand this, let us focus on the simulated trajectories starting around j/q shown 
on figure[6lja) (i? — 0.67, Qi/Qo = 0.2). These trajectories must be analysed in compari- 
son with the unperturbed fiow streamlines, in particular the fiuid separating streamline, 
starting at y — yf and ending up against the front wall at a stagnation point. 

Particles starting around y^ < yf show a clear attraction towards the low flow rate 
branch (displacement along the y axis) as they enter the bifurcation. More precisely, there 
are three types of motions: for low initial position yo (in particular j/o = 0), particles go 
directly into the high flow rate branch. Similarly, above y^, the particles go directly into 
the low flow rate branch. Between some y^* > and yg, the particles first move towards 
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Figure 6. Grey lines: some trajectories of a i? = 0.67 particle when Qi/Qo — 0.2 for (a) 
branches of equal widths, (b) daughter branches 2.5 times wider than the inlet branch and 
(c) daughter branches 7.5 times wider than the inlet branch. The unperturbed fluid separating 
streamline starting at y — yf is shown in black. The particle is shown approximatively at its 
stagnation point. 



the low flow rate branch, but finally enter the high flow rate branch: the initial attraction 
towards the low flow rate branch becomes weaker and the particle eventually follows the 
streamlines entering the high flow rate branch. This non monotonous variation of j/q for a 
particle starting just below j/g is also seen in experiments, as shown in figure|4l right part: 
the third position of the vesicle is characterized by a yo slightly higher than the initial 
one. Back to the simulations, note that, at this level, there is still some net attraction 
towards the low flow rate branch: the particle stagnation point near the opposite wall is 
still below the fluid separating streamline (that is, on the high flow rate side). This two- 
step effect is even more visible when the width 2a of the daughter branches is increased, 
so that the entrance of the bifurcation is far from the opposite wall, as shown on flgures 
[ni[b,c). The second attraction is, in such a situation, more dramatic: for a — 7.5, the 
particle stagnation point is even on the other side of the fluid separating streamline, that 
is, there is some attraction towards the high flow rate branch! Thus, there are clearly 
two antagonistic effects along the trajectory. In the first case of branches of equal widths, 
where the opposite wall is close to the bifurcation entrance, the second attraction towards 
the high flow rate branch coexists with the attraction towards the low flow rate branch 
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and finally only diminishes it. 



These two effects occur in two very different situations. At the entrance of the chan- 
nel, an attraction effect must be understood in terms of streamlines crossing: does a 
pressure difference build up orthogonally to the main flow direction? Near the oppo- 
site wall, the flow is directed towards the branches and being attracted means flowing 
up- or downstream. In both cases, in order to discuss whether some pressure difference 
builds up or not, the main feature is that, in a two-dimensional Stokes flow between 
two parallel walls, the pressure difference between two points along the flow direction 
scales like AP oc Q/h'^, where Q is the flow rate and h the distance between the two 
walls. This scaling is sufficient to discuss in a first order approach the two effects at stake. 



The second effect is the simplest one: indeed, the sphere is placed in a quasi-elongational, 
but asymmetric, flow. As shown on figure [Tl^b), around the fiow stagnation point, the 
particle movement is basically controlled by the pressure difference P2 — Pi, than can be 
written (Pg — P{) — (Pq — P2). Focusing on the y component of the velocity field, which 
becomes all the more important as a is larger than 1, we have Pq — P/ oc Qi/{a — R)^. 
Around the flow stagnation point, the pressure difference P2 — P{ has then the same 
sign as Qi — Q2 and is thus negative, which indicates attraction towards the high flow 
rate branch. For wide daughter branches, when this effect is not screened by the first 
one, this implies that the stagnation point for particles is above the fluid separating line, 
as seen on flg urelBTc). The argument that we use here is similar to the one introduced 
by Fung (see lFung|l993t Yang et al. 2006 ) but resolves only one part of the problem. 
Following these authors, it can also be pointed out that the shear stress on the sphere 
is non zero: in a two-dimensional Poiseuille flow of width h, the shear rate near a wall 
scales as Q/h^, so the net shear stress on the sphere is directed towards the high flow 
rate branch, making the sphere roll along the opposite wall towards this branch. 

Finally, thi s situation is similar to th e one of a flow around an obstacle, that was 
considered in El-Kareh fc Secomb ( 2000( ) as a model situation to understand what hap- 
pens at the bifurcation. Indeed, the authors find that spheres are attracted towards the 
high velocity side of the obstacle. However, we show here that this modeling is mislead- 
ing, as it neglects the first effect, which is the one which eventually governs the net effect. 



This first effect leads to an attraction towards the low flow rate branch. To understand 
this, let us consider a sphere located in the bifurcation with transverse position yo = yf. 
The exact calculation of the flow around it is much too complicated, and simplifications 
are needed. Just as we considered the large a case to understand the second mechanism 
eventually leading to attraction towards the high flow rate branch, let us consider the 
small a limit to understand the first effect: as soon as the ball enters the bifurcation, it 
hits the front wall. On each side, we can write in a first approximation that the flow rate 
between the sphere and the wall scales as Q oc APh^, where AP = Pq — Pi is the pressure 
difference between the back and the front of the sphere, and h the distance between the 
sphere and the wall (see flgure[7^). 

Since the ball touches the front wall, the flow rate Q is either Qi or Q2 and is, by 
deflnition oi yf, the integral of the unperturbed Poiseuille flow velocity between the wall 
and the y = y/ line, so Q (x — /3, where h — 1 ±yf (see flgure[71[a) for notations). 

We have then, on each side: 



A/^oc^^. 



(4.1) 
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Figure 7. (colour online) Scheme of the geometry considered for the two effects occurring in 
the bifurcation, (a) Entrance of the bifurcation: attraction towards the low flow rate branch 
(Pi < P2). (b) opposite wall: attraction towards the high flow rate branch {Pi > P2). 



To make the things clear, let us consider then the extreme case of a flat particle: h = h. 
Then AP oc — 1/3 is a decreasing function of h, that is, a decreasing function of 
Q. Therefore, the pressure drop is more important on the low flow rate side, and finally 
Pi < P2: there is an attraction towards the low flow rate branch. This is exactly the 
opposite result from the simple view claiming that there is some attraction towards the 
high flow rate branch since AP scales as Q/h^ so as Q. Since one has to discuss what 
happens for a sphere in the vicinity of the separating line, Q and h are not indepen- 
dent. This is the key argument. Note finally that there is no need for some obstruction 
arguments to build up a different pressure difference on each side. It only increases the 
effect since the function h (/i^ — h^/3)/{h — R)^ decreases faster than the function 
h i-> {h^ — h^/3)/h^. One can be even more precise and take into account the variations 
in the gap thickness as the fiuid flows between the sphere to calculate the pressure drop 
by lubrication theory. Still, it is found that AP is a decreasing function of h. 

In the more realistic case a ~ 1, the flow repartition becomes more complex, and the 
particle velocity along the x axis is not zero. Yet, as it is reaching a low velocity area 
(the velocity along x axis of the streamline starting at yj drops to 0), its velocity is 
lower than its velocity at the same position in a straight channel. In addition, as the flow 
velocities between the sphere and the opposite wall are low, and since the fluid located 
e. g. between yf and the top wall will eventually enter the top branch by definition, we 
can assume it will mainly flow between the sphere and the top wall. Note this is not true 
in a straight channel: there are no reasons for the fiuid located between one wall and 
the y = yo line, where yo is the sphere lateral position, to enter completely, or to be the 
only fluid to enter, between the wall and the particle. Therefore, we can assume that the 
arguments proposed to explain the attraction towards the low flow rate branch remain 
valid, even though the net effect will be weaker. 

Note finally that, contrary to what discussed for the second effect, the particle rotation 
probably plays a minor role here, as in this geometry the shear stress exerted by the fluid 
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Finally, this separation into two effects can be used to discuss a scenario for bifurcations 
with channels of different widths: if the inlet channel is broadened, the first effect becomes 
less strong while the second one is not modified, which results in a weaker attraction 
towards the low fiow rate branch. If the outlet channels are broadened, as in figures 
[SIb,c), it becomes more subtle. Let us start again by the second effect (migration up- 
or downstream) before the first effect (transverse migration). As seen on figure [51 the 
position of the particle stagnation point (relatively to the flow separating line) is an 
increasing function of a, so the second effect is favoured by the broadening of the outlets: 
for a oo, we end up with the problem of ffow around an obstacle, while for small 
a, one cannot write that the width of the gap between the ball and the wall is just 
a — R, therefore independent from Qi, as it also depends on the y position of the particle 
relatively to yQ. In other words, in such a situation, the second effect is screened by 
the first effect. On the other hand, as a increases, the distance available for transverse 
migration becomes larger, which could favour the first effect, although the slow down of 
the particle at the entrance of the bifurcation becomes less pronounced. 

Finally, it appears to be difficult to predict the consequences of an outlet broadening: 
for instance, in our two-dimensional simulations presented in figure[n](i? = 0.67, Qi/Qo = 
0.2), ?/q varies from 0.27 when the outlet half-width a is equal to 1, to 0.31 when a is 
equal to 2.5 and drops down to 0.22 for a — 7.5! Note that the net effect is always an 
attraction towards the low flow rate branch (j/q < 2//). 

For daughter branches of different widths, it was illustrated in iBarber et al. (|2008h 



that the narrower branch is favoured. This can be explained through the second effect 
(see figurelZla): the pressure drop Pq — PI oc Qi/ (a — R)^ increases when the channel width 
decreases, which favours the narrower branch even in case of equal fiow rates between 
the branches. 



4.2. The consequences on the final distribution 

As there is some attraction towards the low flow rate branch, we could expect some 
enrichment of the low flow rate branch. However, as already discussed, even in the most 
uniform situation, the presence of a free layer near the walls will favour the high flow 
rate branch. We discuss now, through our simulations, the final distribution that results 
from these two antagonistic effects. 

As in most previous papers of the literature, we focus on the case of uniform number 
density of particles in the inlet {n{y) = 1 in equation II. ip . In order to compute the final 
splitting Ni/Nq of the incoming particles as a function of flow rate ratio Qi/Qo one needs 
to know, according to equation (jl.ip . the position j/q of the particle separating line and 
the velocity w* of the particles in the inlet channel. From figure [5] we see that i/q depends 
roughly linearly on {Qi/Qo — 1/2), so we will consider a linear fit of the calculated data 
in order to get values for all Qi/Qq. The longitudinal velocity u* was computed for 
all studied particles as a function of transverse position y^. As shown on figure 13 the 
function wJ(j/o) is well describ ed by a quart i c func tion u%{yQ) = ay^ + /3t/§ -|- 7, which is 
an approximation also used in lBarber et al. ( 20081 ). Values for the fitting parameters for 
this velocity profile and for the linear relationship y^ — £, y- {Qi/Qo — 1/2) are given in 
tabled 

The evolution of Ni/Nq as a function of Qi/Qo for two-dimensional rigid spheres is 
shown on figure |9] for two representative radii. By symmetry, considering Qi < Q2 is 
sufficient. In order to discuss the enrichment in particles in the high fiow rate branch 
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Figure 8. Some longitudinal velocity profiles u*{yo) for two-dimensional spheres of different 
radii R in the inlet channel where a Poiseuille flow of velocity Ux{y) = 1 — is imposed at 
infinity. The full lines show the fits by quartic law uj(yo) ~ aj/o + PVo + 7- 



7? 

a. 
P 
7 



0.25 
-0.96 
-0.85 
0.96 
-1.35 



0.42 
-3.45 
-0.65 
0.91 
-1.25 



0.48 
-4.77 
-0.70 
0.89 
-1.17 



0.53 
-6.33 
-0.64 
0.87 
-1.09 



0.60 
-9.52 
-0.61 
0.84 
-1.01 



0.67 
-11.6 
-0.71 
0.81 
-0.90 



0.71 
-12.6 
-0.73 
0.79 
-0.81 



0.80 



0.75 



Table 1. Values for the fitting parameters (a,/?, 7) for the longitudinal velocity 
"^^(2/0) ~ oiyt + /3j/o 4- 7 of a two-dimensional sphere of radius _R in a Poiseuille fiow of imposed 
velocity at infinity Ux{y) = 1 — y^; for R — 0.80, the velocity profile is too fiat to be reasonably 
fitted by a 3-parameter law, since all velocities are equal to 0.75 ± 0.005 in the explored interval 
yo € [—0.15; 0.15]. We also give the values for the fitting parameter ^ of the linear relationship 
between particle separating line position y^ and flow rate ratios: yo = ^ x {Qi/Qo — 1/2). For 
= 0.8, the strong confinement leads to numerical problems as the sphere approaches the walls. 



(branch 2 then), it is also convenient to consider directly the volume fraction variation 

$2/^0 = (A^2/02)/(iVo/Qo). 

When Qi/Qq — 1/2, the particles flow splits equally into the two branches: A^i — Nq 
and $2 = ^"0- For all explored sizes of spheres, when the flow rate in branch 1 decreases, 
there is an enrichment in particles in branch 2, which is precisely the Zweifach-Fung 
effect: Ni/Nq < Qi/Qq or ^^2/^0 > 1- Then, even in the most homogeneous case, the 
attraction towards the low flow rate branch is not strong enough to counterbalance the 
depletion effect that favours the high flow rate branch. However, this attraction effect 



Spheres in the vicinity of a bifurcation 



21 



0.5 



0.4 - 

0.3 - 

0.2 - 

0.1 - 



•R=0.42 

■ R=0.42 - no-attraction assumption 
-R=0.71 

■ R=0.71 - no-attraction assumption 




Q /Q 

1 

Figure 9. Full lines: spheres relative distribution Ni/No and volume fraction $2/$o as a func- 
tion of flow rate distribution Qi/Qo from our two-dimensional simulations for two representative 
radii R = 0.42 and R = 0.71. The curves are straightforwardly derived from equations (|l.ip and 
1.21 and computed values of y% and u% (table [l}. The results are compared with the hypothesis 
where the particles would follow the streamlines (j/g =2//) (dashed lines). 



cannot be considered as negligible, in particular for large particles: while, in case the 
particles follow their underlying fluid streamline, the maximum emichment in the high 
flow rate branch would be around 40% for R = 0.71, it drops down to less than 17% in 
reality. Similarly, the critical flow rate ratio Qi/Qo below which no particle enters into 
branch 1 is greatly shifted: from around 0.29 to around 0.15 for R = 0.71. For smaller 
spheres (i? = 0.42), this asymmetry in the distribution between the two branches is weak: 
while the maximum enrichment in the high flow rate branch would be around 15% in a 
no-attraction case, it drops to less than 8% due to the attraction towards the low flow 
rate branch. 

When the flow Qi is equal to zero or Qo/2, $2 is equal to $0; thus there is a maximal 
enrichment for some flow rate ratio between and 1/2. The increase in $2 with the 
decrease of Qi (right part of the curves of figure |9]) is mainly due to the decrease of 
the relative importance of the free layer near the wall on the side of branch 2. Two 
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mechanisms are responsible for the decrease of $2 when Qi decreases (left part of the 
curves of figm'e[9]): first, when no particle can enter the low flow rate branch because 
its hypothetical separation line j/q is above its maximum position f — i?, then the high 
flow rate branch receives only additional solvent when Qi/Qo decreases and its particles 
are more diluted. Then all curves fall down on the same curve $2/^0 = 1/(1 ~ Qi/Qo) 
corresponding to A^i = (or N2 = Nq), which results in a sharp variation as Qi/Qo 
goes trough the critical flow rate ratio. A smoother mechanism is also to be taken into 
account here, as it is finally the one that determines the maximum for smaller R. As 
Qi/Qo decreases, branch 2 recruits fluid and particles that are closer and closer to the 
opposite wall. As seen in figure [H the discrepancy between the flow and particle velocities 
increases near the walls, so that A^2 increases less than Q2: the resulting concentration 
in branch 2 finally decreases. 

Finally, for applicative purposes, the consequences of the attraction towards the low 
flow rate branch are twofold: if one wishes to obtain a particle-free fluid (e.g. plasma 
without red blood cells), one has to set Qi low enough so that A^i = 0. Due to attraction 
towards the low flow rate branch, this critical flow rate is decreased and the efficiency 
of the process is lowered. If one prefers to concentrate particles, then one must find the 
maximum of the ^2/^0 curve. This maximum is lowered and shifted by the attraction 
towards the low flow rate branch (see figure[H]). Note that for small spheres (e.g. R — 0.42) 
the position of the maximum does not correspond to the point where Ni vanishes; in 
addition, the shift direction of the maximum position depends on the spheres size: while 
it shifts to lower Qi/Qo values for R = 0.71, it shifts to higher values for R — 0.42. 

The choice of the geometry, within our symmetric frame, can also greatly modify the 
efficiency of a device. Since the depletion effect eventually governs the final distribution, 
narrowing the inlet channel is the first requirement. On the other hand, it also increases 
the attraction towards the low flow rate branch, but one can try to diminish it. As 
discussed in the preceding section, this can be done by widening reasonably the daughter 
branches. For instance, if their half-width is not 1 but 2.5, as in figure iGjb), the slope ^ 
in the law Vq = £, ^ {Qi/Qo — 1/2) increases by around 15% for R — 0.67. The critical 
Qi/Qo below which no particle enters the low flow rate branch increases from 0.13 to 
0.19, which is good for fluid-particle separation, and the maximum enrichment ^2/^0 
that can be reached is 22% instead of 15%. Alternatively, since the attraction is higher in 
two dimensions than in three, we can also infer that considering thicker channels, which 
does not modify the depletion effect, can greatly improve the final result. Note that 
this conclusion would have been completely different in case of high fiow rate branch 
enrichment due to some attraction towards it, as claimed in some papers: in such a 
case, confining as much as possible would have been required, as it increases all kinds of 
cross-streamline drifts. 

4.3. Consistency with the literature 

We now come back to the previous studies already discussed in section [5] in order to 
check the consistency between them and our results. 

The only paper dealing with the position of the particle separating streamline was the 
one bv lBarber et al . (2008), where a symmetric Y-shaped bifurcation is studied (branches 
leaving the bifurcation with a 45° angle relatively to the inlet channel, see figure [TJa)). 
In figure ITUT a) we compare their results with our simulations in a similar geometry. The 
agreement between the two simulations (based on two different methods) is very good, 
except for large particles {R = 0.67) and low Qi/Qq. Note that Barber et al. have chosen 
to consider branches whose widths follow the law = -|- , where wq is the width 
of the inlet branch and wi and W2 are the widths of the daughter branches. This law has 
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Figure 10. (a ) Position of the particle separation line yo in a symmetric Y-shaped bifurcation: 
according to Barb er et al\ (j2008') (data extracted from figure 4 of the cited paper) and accord- 
ing to our simulations. The results for similar spheres in our T geometry are also shown, (b) 
Trajectories from our simulations in the T- and Y- shaped bifurcation, for similar sphere size 
{R = 0.67) and flow rate ratio (Qi/Qo ~ 0.2). Full lines: T geometry; dashed lines: Y geome- 
try. The corresponding separating streamline positions (respectively, yo{T) and VoiY)) are also 
indicated. The sphere that is depicted is located at its stagnation point y'g (see figure |6} in the 
Y geometry. 
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Figure 11. Particle distribution as a function of flow rate ratios for spheres of medium size, 
according to our simulations and according to the simulations shown in .Audet fc Olbricht. (|1987i ') 
(same data as plotted in figure [5]). 



been shown to describe approximately the rela tionship between vessel diameters in the 



arteriolar network (see iMavrovitz fc Rovlll983[ ). With our notations, they thus consider 
a — \/T72 ~ 0.79, while we focused on a = 1 in order to compare with the T-shaped 
bifurcation. In addition, their apex has a radius 0.75 (for the R = 0.67 case) while ours 
is sharper (radius of 0.1). These differences seem to impact only partly the results, as 
discussed above. We can expect this slight discrepancy to be due to the treatment of the 
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numerical singularities that appear when the particle is close to one wall. For R — 0.67, 
the maximum position ?/o is 0.33, which is close to the separating streamline position. 

It is also interesting to compare our results in the Y-shaped bifurcation with the results 
in the T geometry, which was chosen to make the discussion easier. We can see that, for 
low enough Qi/Qo, the attraction towards the low flow rate branch is slightly higher. 
This can be understood by considering a particle with initial position tjq slightly below 
the critical position yQ found in the T geometry: in the latter geometry, it will eventually 
enter the high flow rate branch, by definition of j/q . As shown in figure [TUT b) , in the Y 
geometry, this movement is hindered by the apex since the final attraction towards the 
high flow rate branch occurs near the opposite wall (the second effect discussed in section 
I4.1.2[) . Finally from this comparison we see that comparing results in T and symmetric 
Y geometry is relevant but for highly asymmetric flow distributions. 



In sec tion [2l the analysis of t he two-dimensional simulations for R — 0.5 spheres 
shown in lAudet fc Olbrichtl (jl987l) showed that there should be some attraction towards 
the low fiow rate branch. Our simulations for R = 0.48 showed that this effect is non 
negligible (figure [SId) and modifies greatly the final distribution (figure |9|). Finally, we 
can see in figure [TT] that our simulations give similar results as the simulation by Audet 
and Olbricht. 

As for the experiments presented in Yang et al. ( 2006[) for R — 0.46, we showed that 
the final distribution was consistent with a no-attraction assumption. As we showed in 
figure [5ja), in a three-dimensional case, the attraction towards the low fiow rate region 
is weak for spheres of radius R ~ 0.5 or smaller, which is again coherent with the results 
of Yang et al.. Note that, while their results were considered by the authors as a basis to 
discuss some attraction effect towards the high fiow rate branch, we see that their final 
distributions are just reminiscences of the depletion effect in the inlet channel. 

The other consistent set of studies in the literature deals with large balls in three 
dimensional channels. We have studied balls of radius R = 0.71 that stop entering branch 
1 when Qi/Qo ^ 0.22 (figure[5^), while this critical flow rate would be around 0.29 in case 
they would follow the fluid streamlines. This critical flow rate is expected to be slightly 
higher for larger balls of radius R ~ 0.8 , but far lower than 0.35, which would be the 
no-attraction case. In the experiments of lRoberts k, Olbricht ( 2006 ). some balls are still 
observed in branch I when Qi/Qo ~ 0.22 (figure[2l), indicating a stronger attraction effect 
towards the low fiow rate branch, which can be associated to the fact that the authors 
considered a square cross section channel, while the confinement in the third direction 
is 0.5 < 0.71 in our cas e. The experiment s wit h circular cross section chan nels lead to 
contradictory results: in lChien et al. ( 1985 ) and Ditchfield Sz Olbricht (1996), the results 
were consistent with a no-attracti on assumption, there f ore th ey are in contradiction with 
our results. On the contrary, in Roberts fc Olbricht ( 20031) . the critical flow rate for 
R — 0.77 is around 0.2, which would show a stronger attraction than in our case. Note 
that all these apparently contradictory observations are to be considered keeping in mind 
that the data of Ni/Nq as a function of Qi/Qo are sometimes very noisy in the cited 
papers. 



5. Conclusion 

In this paper, we have focused explicitly on the existence and direction of some cross 
streamline drift of particles in the vicinity of a bifurcation with different flow rates in the 
daughter branches. A new analysis of some previous unexploited results of the literature 
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first gave us some indications on tlie possibility of an attraction towards tlie low flow 
rate branch. 

Then the first direct experimental proof of attraction towards the low flow rate branch 
was shown and arguments for this attraction were given with the help of two-dimensional 
simulations. In particular, we showed that this attraction is the result of two antagonistic 
effects: the first one, that takes place at the entrance of the bifurcation, induces migration 
towards the low flow rate branch, while the second one takes place near the stagnation 
point and induces migration towards the high flow rate branch but is not strong enough, 
in standard configurations of branches of comparable sizes, to counterbalance the first 
effect. 

This second effect is the only one that was previously considered in most papers of the 
literature, which has lead to the misleading idea that the enrichment in particles in the 
high flow rate branch is due to some attraction towards it. On the contrary, it had been 
argued by Barber et al. that there should be some attraction towards the low flow rate 
branch. By distinguishing the two effects mentioned above, we have tried to clarify their 
statements. 

In a second step, we have discussed the consequences of such an attraction on the 
final distribution of particles. It appears that the attraction is not strong enough, even 
in a two-dimensional system where it is stronger, to counterbalance the impact of the 
depletion effect. Even in the most homogeneous case where the particles are equally 
distributed across the channel but cannot approach the wall closer than their radius, the 
existence of a free layer near the walls favours the high flow rate branch, which eventually 
receives more particles than fluid. 

However, these two antagonistic phenomena are of comparable importance, and none 
can be neglected: the particle volume fraction increase in the high flow rate branch is 
typically divided by two because of the attraction effect. On the other hand, the initial 
distribution is a key parameter for the prediction of the final splitting. For deformable 
particles, initial lateral migration can induce a narrowing of the ir distribution , which will 
eventually favours the high flow rate branch. For instance, in iBarber et al\ 12008), the 
authors had to adjust the free layer width in their simulations in order to fit experimental 
data on blood fiow. On the other hand, in a network of bifurcations, the initially centered 
particles will find themselves close to one wall after the first bifurcation, which can favour 
a low flow rate branch in a second bifurcation. 



Note finally that, as seen in lEnden fc Popell (|l992l ). these effects become weaker when 



the confinement decreases. Typically, as soon as the sphere diameter is less than half the 
channel width, the variations of volume fraction do not exceed a few percent. 

For applicative purposes, the consequences of this attraction have been discussed and 
some prescriptions have been proposed. Of course, one can go further than our symmetric 
case and modify the angle between the branches, or consider many-branch bifurcations, 
and so on. However, the T-bifurcation case allowed to distinguish between two goals: 
concentrating a population of particles, or obtaining a particle-free fluid. The optimal 
configuration can be different according to the chosen goal. Similar considerations are 
also valid when it is ab out doing som e sorting in polydisperse suspensions, which is an 



important activity (see iPammd 120071 ): getting an optimally concentrated suspension of 
big particles might not be compatible with getting a suspension of small particles free of 
big particles. 



Now that the case of spherical particles in a symmetric bifurcation has been studied 
and the framework well established, we believe that quantitative discussions could be 
made in the future about the other parameters that we put aside here. In particular, 
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discussing the effect of the deformability of the particles is a challenging problem if 
one only considers the final distribution data, as the deformability modifies the initial 
distribution, but most probably also the attraction effect. In a network, the importance 
of these contributions will be different according, in particular, to the distance between 
two bifurcations, so they must be discussed separately. 

Considering concentrated suspensions is of course the next challenging issue. Particles 
close to each other will obviously hydrodynamicaly interact, but so will distant particles, 
through the modification of the effective resistance to fiow of the branches. In such a 
situation, considering pressure driven or flow rate driven fluids will be different. 

For concentrated suspensions of deformable particles in a network, like blood in the 
circulatory system, the relevance of a particle-based approach can be questioned. Histori- 
cal models for the major blood flow phenomena are continuum models with some ad-hoc 
parameters, which must be somehow re lated to the i ntrins ic mechanical properties of 
the blood cells (for a recent example, see Obrist et al. ( 2010l )). Building up a bottom- up 
approach in such a system is a long quest. For dilute suspensions, some links between the 
microscopic dyn amics of lipid vesicles and the rheology of a suspension have been rec ently 
estabhshed (see banker fc Misbahll2007t IVitkova eraZll2008t IChigliotti eFaZlbOlOf). For 
red blood cells, that exhibit qualitatively similar dynamics (see Abkarian et al\ 120071 : 
Deschamps et al. 20091 : Noguchi 2010l : Farutin et al. 20 l O': Dupirc et aZ."2010'), we can 
hope that such a link will soon be established, following iVitkova et al.. (.20083 . For con- 
fined and concentrated suspensions, the distribution is known to be non homogeneous, 
which has direct consequences on the rheology (the Fahraeus-Lindquist effect). Once 
again, while empirical macroscopic models are able to describe this reality, establishing 
the link between the viscosity of the suspension and the local dynamics is still a challeng- 
ing issue. The final distribution of the fiowing bodies is the product of a balance between 
migration towards the center, which has already been discussed in the introduction of 
the present paper, and interactions betwee n them that can broaden the distribution (see 
Kantsler et al. [ boost IPodgorski et a/.|[2O10[) . The presence of d e forma ble boundaries also 
needs to be taken into account, as shown in iBeaucourt et al\ (j2004l ). In the meantime, 
the development of simulations techniques for quantitative three- dimensional approaches 
is a crucial task, w hich is becoming more and more feasible (see McWhirter et al. 20091 : 
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Appendix. 

In this appendix, details for the derivation of equations p.l3p - p.l4p - (l3.15p from equa- 
tions (pm| - ((XTT|) - (IXT^ are given. 

We introduce first the scalar product in Lp'ip.fY as follows: 

Vf,geL2(%)2, <f,g>i2(o,)2= / f.g. 

The variational formulation of problem (I3.10p - (l3.1ip - (l3.12p is obtained by taking the 
scalar product of the equation p.lOp in L^{Q,f Y with a test function v G i?Q(r2/)^ and 
we multiply equation p. lip by a test function q £ io(^)- It leads to this problem: find 
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(u,p) e H^iflf)^ X Lliflf) such that: 

-2i^ [ (V-t(u)).v+ / Vp-v = 0,Vvei/i(%)2, (Al) 
JUf Jet] 

f gV-u = 0,VgeL2(%), (A 2) 

u = fon9%. (A3) 
Applying Green's formula to equation (jA ip we obtain 

2v t(u) : Vv - 2i^ / T(u)n • v 

Jn, JdUf 



pV • V + / pv • n = 0. (A 4) 

Where n denotes the outer unit normal on dflf. Taking into account that v vanishes 
on 9ri/ (recall that we have chosen the test function v G iJQ(51/)^), the problem (jA ip - 
(|A2p - (|A3p is now equivalent to this one: find {u,p) e H^{^fY x -^o(^/) ^^^^ that: 

2v ( t(u):Vv-/" pV-v = 0,Vvei7i(%)2, (A 5) 

[ qV ■u = 0,VqeLl{nf), (A 6) 

u = fondnf. (A 7) 

Note that t(u) is symmetric (t(u) : Vv = t(u) : (Vv)*). So that we can write t(u) : 
Vv = r(u) : r(v). Finally, the variational formulation of our initial problem p.ll) - p.2p - 
(|33)) is given by: find (u,p) e H'^i^ff x Lg(rJ/) such that: 

"(%)^ (A8) 



2v ( t(u):t(v)-/" pV-v = 0,Vvei7o^( 

[ qV ■u = o,yqeLl{nf), 



(A 9) 



u = fonarj/. (A 10) 

Remark 2. we have 

r(u) : t(v) = t(u) : Vv = i Vu : Vv + ^(Vu)* : Vv, (A 11) 



the first integral in equation HA ($)) can he rewritten thanks to this identity 

I T(u):T(v) = i/ Vu:Vv. (A 12) 

Jnf 2 

Indeed, by integration by part and using the incompressibility constraint ^ • u — we 
tiave / (Vu)* : Vv = 0. Thus we can retrieve the formulation of our problem as a 

minimization of a kind of energy. The velocity field u is then the solution of this problem 

J(u) ^ inf J(v), (A 13) 
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where 

J(v) = / Vv : Vv = /" r(v) : t(v). (A 14) 
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